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A hybrid model of molecular dynamics and continuum mechanics is introduced to study a system 
of vertically shaken granular layers. Despite the simplicity the model shows pattern formation in 
the granular layers due to the formation of heaplets. We show from a simple analysis that the onset 
of pattern formation is density dependent and this result is justified by the subsequent computer 
simulations. Our simulations also show that the heaping process can be divided into three stages: 
an early stationary stage, an intermediate growing stage and a late-time saturated stage. In the 
early stage, the average volume of the heaplets remains almost unchanged until the system crosses 
over to the intermediate growing stage. The length of time that system remains in the early stage 
defines the onset time of the instabilities, ko which depends on the shaking intensity, 7. In the 
growing stage, the average volume of the heaplets grows with time and can be approximated by a 
power law with a shaking intensity dependent growth exponent, z. From our simulation, we show 
that the growth exponent, z depends linearly on the shaking intensity, 7. There is a critical shaking 
intensity -y c ~ 14, such that for shaking intensities greater than j c , z cannot be properly defined. 
The late-time saturated stage is where most of the particles are trapped in a big heap and this big 
heap is in equilibrium with the surrounding granular gas. 

PACS numbers: 61.43.Gt,45.70.Qj,83.80.Fg 
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I. INTRODUCTION 

We are interested in pattern formation in a thin layer of 
granular particles upon a plate. We begin by distinguish- 
ing between experiments on systems in which the plate is 
vertically vibrated |jj or in which the plate is tapped ||] 
or shaken. In the vertically vibrating systems, the par- 
ticles dance upon the oscillating plate and the surface of 
the layer becomes patterned with regions characterised 
by different modes of oscillation or by the the same modes 
of oscillation but different phases. The features described 
as oscillons || and the interface between regions Q are 
both of interest here. In the tapped or shaken system 
there is a long period between taps/shakes in which the 
plate is at rest and the particle layer relaxes in this pe- 
riod to a static configuration in which the surface of the 
layer is marked by ripples or heaplets. In other words 
the perturbations are discrete and are followed by peri- 
ods of stasis in distinction to the vibrated system which 
is continuously moving. 

Following the original observation of square and 
striped patterns by Melo et. al |) in the vibrat- 
ing system, there have been many experimental observa- 
tions H ||, and theoretical discussions j|, ||, || 
for this type of pattern formation. All these works and 
many others have contributed greatly and also establish 
a good understanding of the dynamics of the pattern 
formation of granular layer under continuous vibration. 
On the other hand in the tapped/shaken system, Duran 
has observed the formation of ripples JL^ | and granular 
heaplets Q in a tapped granular layer. However, there is 
less literature available on this class of experiments, and 
our paper is an attempt to fill this gap. 

We now further distinguish between vertically shaken 



granular layers and tapped granular layers. In the shaken 
system the plate moves up and then is withdrawn sharply 
to a neutral position. When the plate reaches the high- 
est point, the particles continue their motion and then 
travel ballistically, falling back to the stationary plate, 
rebounding and then settling down and relaxing into a 
static configuration. The plate is subsequently moved up 
again and the process repeated. The lateral redistribu- 
tion, which produces the patterns, takes place when the 
falling particles rebound from the plate and relax. On 
the other hand in the tapped system the plate scarcely 
moves but, as in the Newton's cradle experiment, the 
tap is imparted impulsively to the particles which jerk 
off the plate, fall back, rebound and relax. Most of the 
lateral redistribution takes place at the jerk. However, 
shaken and tapped granular layers also share some com- 
mon properties such as the formation of isolated granular 
heaplets, which coarsen with the number of taps/shakes, 
and the rate of coarsening of the heaplets increases as 
the intensity of the shocks increases. We have discussed 
tapped layers in a previous article |l3| |. The present is 
concerned mainly with shaken granular layers. 

This article is organised as follows: first, we introduce a 
simple model to study shaken granular layers. The model 
is based on a hybrid of a molecular dynamics description 
and continuum modelling of a granular sheet. Then we 
show typical patterns produced by the model and these 
patterns resemble typical experimental morphologies of 
shaken granular layers with isolated heaplets. An analy- 
sis is performed and leads to the prediction that the de- 
pinning shaking intensity of the granular layers is density 
dependent, which is later justified from the simulations. 
Later we show from the simulation that the coarsening 
of the pattern comes in three stages: an early stationary 
stage, a growing stage and a saturated stage. In the early 
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stage, the layer remains fiat, with small fluctuation su- 
perimposed, for a short duration kg which depends on the 
shaking intensity 7. The early stage is characterised by a 
small value V (defined later), which measure the degree 
of heaping. Then we discuss the cross-over time from 
the early stage to the growing stage, which indicates the 
onset time ko of instabilities. In the growing stage, we ap- 
proximate the growth of the average volume of heaplets 
by a power law. The growth exponent z is 7-dependent 
and the power law approximation fails when 7 is larger 
than a critical value, j c . Finally, the late-time saturated 
stage is where all the granular heaplets are merged into a 
big heap and the big heap is in equilibrium with a scat- 
tering of particles on the plate which can be described as 
a granular gas. The time of onset of saturation is depen- 
dent on the size of the system. In our simulations, the 
system size is small (with periodic boundary conditions), 
and as a result the time of onset will almost certainly be 
small compared to that experimentally observed in real 
systems. 



II. MODEL 

The spatial distribution of a granular layer can be 
described by a continuum field variable n(r, t), where 
r = {x,y) are the Cartesian coordinates parallel to the 
plate. The number density n(f, t) also represents the 
layer thickness if we take the average diameter of grains 
as unity assuming no compactification throughout the 
whole shaking process [l3|] . When the system is shaken 
the layer moves off the plate to height h(r, t) measured 
from the stationary neutral position of the plate. There is 
also an associated velocity variable u(f,t) = dh(f,t)/dt 
and u(r, t) then represents the average vertical velocity of 
the grains at the point r. As in Rothman jl4| the velocity 
field at each point is coupled dissipatively with the field 
at neighbouring points through a Laplacian, and the dy- 
namics can be described by a generalisation of Newton's 
second law, 



du 
~dt 



vV 2 u — g + B(h, u, a). 



(1) 



Here v is a dissipative coefficient; it represents the lateral 
exchange of velocities when particles collide with each 
other and g is the acceleration due to gravity. The func- 
tion B specifies how the layer bounces when it hits the 
plate. The bouncing term B is zero except when h = 
and it is defined by the mapping 

u — > ~au, (2) 

where a is the coefficient of restitution. The same equa- 
tions here have been used by Rothman [Q and Moon et 
al. |ll| to describe the phenomena of oscillons and stripes 
structures in a vertically vibrated sand layer. 

The coefficient of restitution a is expected to depend 
on n because the energy dissipation is greater in regions 



of higher local density. However, the precise form of a(n) 
is not known, except that it should be a monotonic de- 
creasing function of n. In fact, which we will show in 
the following section, it is possible to deduce roughly the 
functional form of a(n) from experiments by measuring 
the depinning shaking intensities for homogeneous gran- 
ular layer with different thickness. In this paper, for the 
sake of definiteness, a(n) is chosen to be a Lorentzian, 
i.e. 



a{n) 



(3) 



where no and n\ are parameters. 

In the shake the plate moves up and is then withdrawn 
sharply. The particles continue in free flight, all with the 
same velocity, and fall under gravity until they hit the 
plate at rest, all with the same velocity —U. They then 
rebound with local velocity u(f,t) ~ a(n)U which de- 
pends on the local density n(f, t) because of the density 
dependence of the coefficient of restitution a. It is use- 
ful to use 7 = U/Uq as a dimensionless measure of the 
shaking intensity, here Uq — 3cm/s being a very crude 
estimate of a typical take-off velocity of the grains when 
the plate is withdrawn. The purpose of Uq is to allow us 
to present the data on a reasonable scale. 

After the shake, the particles undergo motion as given 
by Eq. ([j]). When the particles at a point fall back to 
the plate again, they bounce and redistribute onto the 
neighbourhood. To describe this, we discretise the sys- 
tem into L x L squares of side a with the square at fi 
containing Ni — n(fi, t) a 2 particles. Then we implement 
the redistribution in the following way: 

1. We define an integer Ni equal to the the next in- 
teger greater than n(fi,t), which is the number of 
particles on the rebounding site. 

2. We surround the rebounding site fi by a disc of 
radius R ma x pi , which is given by 



Rjnax £ 



u 



(4) 



where £ is a parameter which controls the distance 
travelled by the particles. 

3. We specify a hopping angle <j) by choosing a random 
number in the range [0, 2tt] with equal probability. 

4. We specify a hopping range R by choosing a ran- 
dom number in the range [0, R ma x] with equal prob- 
ability. 

5. The angle <f> and hopping range R together define 
a receiving site r'j such that \r'j — fi\ = R and 
the angle between the vector r' j — f and a fixed 
reference direction is d>. 



6. The process is repeated Ni times but receiving sites 
may be repeated and the rebounding site itself may 
be a receiving site. 
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7. On each repetition (except the last) the density 
n{fi,t) on the rebounding site is reduced by unity 
and the density on the receiving site is increased 
by unity. On the last repetition, leftovers are redis- 
tributed. 

The redistribution rules can be summarised by first 
defining 2Vj and the size of the leftovers, i.e. 
Ant = n(f*i,t) — Ni + 1, then repeat the following JV, 
times. 



n(ri,t) 
n(r'j,t) 

where 



n{?i,i) - Aj 



f 

An,- 



1,2,. 
Ni. 



fori = 1,2,. 



,N, 
(5) 

(6) 



The process conserves the number of particles in the re- 
distribution as it should. If R m ax is greater than system 
size the physical reality is not well-represented. 

The mass transfer after redistribution can be describe 
by a continuity equation dtn = — V • J. The constitutive 
current equation is given by 



-rj({3\Vn\ 2 - l)Vn, 
0, 



|Vn| > ^ 



(7) 



where 1/a//5 = tan# c is the critical slope of the heaplet 
side and r] sets the diffusion rate. Here the angle of re- 
pose 9 C is chosen to be 7r/6. The constitutive equation 
is constructed in such a way that when the gradient is 
greater than a critical slope, mass current will flow down 
hill to smooth out density fluctuation, but no mass is 
transfered if the gradient is less than the critical slope, 
because this corresponds to a stable situation. 



IV. RESULTS 



A. Morphology 



Fig. H shows typical results of the pattern formation 
in two different simulations. The top panel corresponds 
to a lower shaking intensity (7 = 4.0), and the lower 
panel corresponds to a higher shaking intensity (7 = 7.0). 
Bright regions correspond to a higher density region of 
the granular heaplets in the granular layer. Dark regions 
correspond to regions of low densities of freely moving in- 
terstitial granular particles (granular gases) . Both panels 
show an effective surface tension as suggested in refer- 
ence [^3). The effective surface tension depends on the 
shaking intensity. In the lower panel a greater effective 
surface tensions is associated with a higher shaking in- 
tensity, and this greater surface tension produces more 
compact heaplets, which in agreement with the findings 
in reference [fill . 




FIG. 1: Two simulations of a bouncing granular layer corre- 
sponding to different tapping intensities. The top panel cor- 
respond to a lower shaking intensity 7 = 4.0, and the lower 
panel corresponds to a higher shaking intensity 7 = 7.0. Each 
column is taken at a specific time and from left to right the 
corresponding number of shakes is k = 5, 10, 15, 30 and 50. 



III. SIMULATION 



The model is studied on a L x L square lattice with pe- 
riodic boundary condition. Initially the granular layer is 
prepared with the density n randomly distributed about 
an average value n. At the beginning of each shake, 
each site is assigned a velocity u(r,t) — a(n{r,t))"]U . 
We choose the parameters hq and n\ in Eq.(|) for a(n) 
such that n,Q = 1/3 and n\ = 2/3. The value of the 
lattice spacing a is set equal to the average diameter 
of the particles and the hopping parameter £ is set to 
7b\fiUn I 2q Jl3fl . The reason we choose £ differently from 
reference |13| and Uq = 3 cm/s is to allow us present data 
on a reasonable scale. Then Eq. (jl|) and the equation of 
continuity with Eq.(Q) are solved using a Runge-Kutta 
method. On each shaking cycle, these equations are it- 
erated until no more lateral mass transfer occurs at each 
site because of the rebound and because of the local gra- 
dient on each site is less than the tangent of the angle of 
repose, 8 C = tt/6 with 1% cutoff. 



B. Depinning Shaking Intensity 

In order to form heaplets in the granular layer, the 
shake intensity not only needs to be strong enough to lift 
grains off the plate but it must also be able to redistribute 
grains so the hopping range must greater than the lattice 
spacing a. One can predict the shaking intensity 70 for 
such a depinning process from a mean field argument. 
Let the thickness of the granular layer be homogeneous 
throughout the system and be equal to n. After the first 
shake, the rebound velocity of the layer is u = a(n)jUo 
everywhere. We can obtain the depinning shaking inten- 
sity 70 from Eq.(||) and setting R m ax(n) — a, so that 

lQ = C(n 2 +nl) (8) 
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and 




However, one might expect the relation given by Eq.(|^) 
to be not strictly quadratic in h when fluctuations are 
important. Fig. @ shows the variation of 70 against aver- 
age density n from the simulation. The graph show that 
the relation between 70 and n is in good agreement with 
the prediction. From the inset of Fig. || shows that the 
exponent is only slightly less than 2 (~ 1.91). This result 
is due to the specific form for the effective coefficient of 
restitution a(n) we chose in Eq.(|J). However, this also 
suggests a possibility that one may construct the func- 
tional form of a(n) from experiments. This can be done 
by noting the relation between depinning shaking den- 
sity 70, which can be measured from experiments, and 
the corresponding average density of layer n is given by 
70 = (a/0 1 ^ 2 [ a (^)]~ ■ We present this as a challenge 
to experimentalists. 




(i - 



I 1 I 1 I 1 I 1 I 1 L 

1 2 3 4 5 



heaplet volume V for different values of 7 as a function 
of the number of shakes, k. 

In each plot in Fig. ||| there are three different stages 
and each stage corresponds to a different underlying pro- 
cess. In the early stage, the volume remains almost un- 
changed in time. In this stage, the fluctuations in the ran- 
domly distributed granular layer are not strong enough to 
trigger coalescence. As a result, the layer remains nearly 
flat and has a small value of V. The second stage is the 
heaplet formation stage, which corresponds to the grow- 
ing regions in the plots. In this stage, high density regions 
of the granular layer act as a kinetic energy sinks, many 
random hopping granular particles hop into these sinks 
and are trapped there. As a result of capturing particles, 
these high density regions grow and merge to form big- 
ger heaplets. Finally, the late time stage is where all the 
the heaplets have joined into a single big heap and stop 
growing. There are some remaining individual particles 
(granular gas) hopping randomly about and occasionally 
encountering the big heap and getting trapped there. Oc- 
casionally particles on the big heap will leave the heap 
and join the surrounding granular gas. When the num- 
ber of particles leaving the big heap is balanced by those 
entering the big heap, the system is now in an equilib- 
rium situation which very similar to that of droplets in a 
supersaturated vapour. However there are very few indi- 
vidual particles involved in the motion across boundaries 
during equilibrium and therefore the fluctuations of the 
volume of big heap are small as observed in Fig. || and 
Fig. ||. Of course the size of the saturated heap depends 
on the size of the simulated system and the average num- 
ber density of particles. Larger system sizes will result 
in a later saturation and bigger saturated heaps. 



FIG. 2: Variation of 70 for different values of n. The inset 
shows the log-log plot for 70 against n. The linear regression 
fit line for the inset gives a slope of 1.91. 



C. Average volume of heaplets 

We are also interested in how the size of heaplets grows 
in time. The relevant length scale can be obtained as the 
average of the thickness at site i weighted by its local 
volume. We know that the thickness is rij and the local 
volume clement of site i is V7 = m<J, where site area a is 
the square of the lattice constant, then 



/ 



Ei hi Vi _ Ei' 



(9) 



Hence the volume of the heaplet is V = I 3 and Fig. ^ 
and Fig. || shows the typical linear and log-log plots of 
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FIG. 3: Log-Log plot of heaplet size, V against number of 
shakes, k for different values of 7. From left to right the 
corresponding value of 7 for each curve is 6, 8, 10, 12 and 
14. Each curve is averaged over 10 different runs. Note the 
break-down of power-law growth in the intermediate region 
for the curves with 7 > 14. 
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D. Growth Exponents 

In the growing stage, the growth exponent z is 7- 
dependent. Fig.([|) shows the growth exponent z for 
different values of 7. As we can see there is a trend 
of increasing growth exponent as the shaking intensity 
7 is increased. The relation between the growth expo- 
nent and the shaking intensity is roughly linear and is 
given by z ~ 0.3487 — 0.644, as shown as the straight 
line in Fig.(^). However, the relation breaks down when 
7 > 14. The reason of the break-down of the power law 
approximation will be discussed in the following section. 
We can understand the dependence of z to 7 from the 
fact that larger shaking intensities give larger kinetic en- 
ergy to allow the system explore greater regions of phase 
space. As a result, systems with larger shaking inten- 
sities are able to find stable configurations easier than 
systems with lower shaking intensities. Of course, the 
final stable configuration is a single heap which favour 
greater dissipation. 



mation fails for 7 > j c . This can be clearly seen from 
Fig.^. Each panel of Fig. ^ contains time series of V for 
10 different runs where each run has a different initial 
configuration but the same shaking intensity. The top 
I 
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FIG. 5: Variation of log 10 fco against log 10 7. It is obvious 
that fco increases with 7. 




FIG. 4: The growth exponent z as a function of shaking in- 
tensity 7. The straight line is the linear regression fit for the 
data, the slope of straight line is 0.348. 



E. Onset Time 

We also note from Fig. ||| that the onset of the growing 
stage occurs at different times for different 7's. Generally 
the onset time fco increases as 7 increases. Fig.(||) plots 
fco against 7 and shows that fco is an increasing function 
of 7. However, for 7 greater than a critical value 7 > 14, 
heaplets can start forming at any time hence no proper 
fco can be defined. The same reason that fco loses its 
predictability also explains why the power law approxi- 



is clear that there is a well defined onset time of instabil- 
ities at fc ~ 6 and consistent growth after the instability 
starts. The irregular slow-down of the growth during the 
cross-over from growing stage to saturated stage is due 
to the formation of multiple metastable heaps and the 
merging process of these heaps is slow. On the other 
hand, curves in the lower panel do not have a well de- 
fined onset time, instabilities in each curve appear at a 
different time. In fact, some curves remain trapped in 
the metastable early stage and instabilities never appear 
through out the simulations. 



V. CONCLUSION 

We have introduced a simple model to study vertical 
shaken granular layers. The model shows the shaking in- 
tensity for the depinning of granular layer is related to 
the layer thickness and to the effective coefficient of resti- 
tution of the granular layer. The dynamics of heaplet 
formation in shaken granular layer can be divided into 
three stages: an early stationary stage (7-independcnt), 
an intermediate growing stage (7-dependent), and a sat- 
urated/equilibrium stage. In the growing stage, the 
growth exponent z is an increasing function of 7. The 
onset time fco increases as 7 increases up to a critical 
value of 7c (— 14) where heaplet growth can start at any 
time dependent on details of the initial configuration. 
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